<?php
/**
	Update History:
	Ver.    TRB#        Date          Author    Note
	1                   2010-09-15    Trungdt   tạo mới
	**/
/*
* Tài liệu tham khảo: http://www.alglib.net/interpolation/polynomial.php
* Module này sử dụng lại gói JAMA 
* @package JAMA
* @author Michael Bommarito
* @author Paul Meagher
* @version 0.1
*
* Function to fit an order n polynomial function through
* a series of x-y data points using least squares.
*
* @param $X mảng các phần tử X
* @param $Y mảng các phần tử Y
* @param $n int bậc tối đa của đa thức. Nên để $n=3 (bậc 3)
* @returns một object với cấu trúc như sau
*    ret_object
*       -> m: kích thước của mảng các hệ số được trả về
*             nếu tham số bậc đầu vào là 3, m sẽ là 4 (vì có bậc 0)
*       -> A: mảng chứa hệ số
*           ->A[$i][0]: giá trị hệ số
*           Chứa các hệ số từ thấp đến cao
*           A[0]: bậc 0, A[1]: bậc 1; A[2]: bậc 2
*/
require_once "Matrix.php";
function polyfit($X, $Y, $n) {  
  for ($i = 0; $i < sizeof($X); $i++) 
	  for ($j = 0; $j <= $n; $j++) 
		$A[$i][$j] = pow($X[$i], $j);
  for ($i=0; $i < sizeof($Y); $i++)
	$B[$i] = array($Y[$i]);   
  $matrixA = new Matrix($A);
  $matrixB = new Matrix($B);
  $C = $matrixA->solve($matrixB);
  return $C->getMatrix(0, $n, 0, 1);
}

function printpoly( $C = null ) {
  for($i = $C->m - 1; $i >= 0; $i-- ) {
	$r = $C->get($i, 0);
	if ( abs($r) <= pow(10, -9) )
	  $r = 0;
	if ($i == $C->m - 1)
	  echo $r . "*B3^{$i}";
	else if ($i < $C->m - 1)
	  echo " + " . $r . "*B3^{$i}";
	else if ($i == 0)
	  echo " + " . $r;
  }
}
//$X=array(301,303,307,309,311,319,321,324,326,329,330,331,332,333,335,337,339,340,343,344,348,350,351,354,355,356,357,361,362,363,365,368,370,371,372,373,375,376,379,380,381,386,387,388,389,390,391,392,393,394,398,399,400,401,404,405,406,408,410,411,412,413,414,417,418,420,421,422,423,424,425,426,427,430,431,432,433,434,435,436,437,438,439,440,441,442,444,445,446,447,448,449,450,451,452,454,456,458,460,461,462,463,464,465,466,467,468,469,470,474,475,477,478,479,481,483,484,486,488,489,493,504,505,507,508,509,512,513,514,515,517,518,519,520,521,522,523,524,525,526,527,528,531,532,533,534,535,536,537,538,544,546,548,550,551,556,558,559,561,565,576,579,585,587,588,590);
//$Y=array(110.1425698,109.2953193,107.6008183,106.7535678,105.9063173,102.5173153,101.6700648,100.399189,99.5519385,98.28106275,97.8574375,97.43381225,97.010187,96.58656175,95.73931125,94.89206075,94.04481025,93.621185,92.35030925,91.926684,90.232183,89.3849325,88.96130725,87.6904315,87.26680625,86.843181,86.41955575,84.72505475,84.3014295,83.87780425,83.03055375,81.759678,80.9124275,80.48880225,80.065177,79.64155175,78.79430125,78.370676,77.09980025,76.676175,76.25254975,74.1344235,73.71079825,73.287173,72.86354775,72.4399225,72.01629725,71.592672,71.16904675,70.7454215,69.0509205,68.62729525,68.20367,67.78004475,66.509169,66.08554375,65.6619185,64.814668,63.9674175,63.54379225,63.120167,62.69654175,62.2729165,61.00204075,60.5784155,59.731165,59.30753975,58.8839145,58.46028925,58.036664,57.61303875,57.1894135,56.76578825,55.4949125,55.07128725,54.647662,54.22403675,53.8004115,53.37678625,52.953161,52.52953575,52.1059105,51.68228525,51.25866,50.83503475,50.4114095,49.564159,49.14053375,48.7169085,48.29328325,47.869658,47.44603275,47.0224075,46.59878225,46.175157,45.3279065,44.480656,43.6334055,42.786155,42.36252975,41.9389045,41.51527925,41.091654,40.66802875,40.2444035,39.82077825,39.397153,38.97352775,38.5499025,36.8554015,36.43177625,35.58452575,35.1609005,34.73727525,33.89002475,33.04277425,32.619149,31.7718985,30.924648,30.50102275,28.80652175,24.146644,23.72301875,22.87576825,22.452143,22.02851775,20.757642,20.33401675,19.9103915,19.48676625,18.63951575,18.2158905,17.79226525,17.36864,16.94501475,16.5213895,16.09776425,15.674139,15.25051375,14.8268885,14.40326325,13.979638,12.70876225,12.285137,11.86151175,11.4378865,11.01426125,10.590636,10.16701075,9.7433855,7.201634,6.3543835,5.507133,4.6598825,4.23625725,2.118131,1.2708805,0.84725525,4.75E-06,-1.69449625,-6.354374,-7.62524975,-10.16700125,-11.01425175,-11.437877,-12.2851275);
//$X = array(1.02, 1.33 , 1.58, 1.78 , 1.97, 2.16, 2.35, 2.52, 2.71, 2.88, 3.07, 3.21, 3.42, 3.63, 3.89, 4.16);
//$Y = array(0,25,50,75,100,125,150,175,200,225,250,275,300,325,350, 375);
$X=array(1,1.5 ,2 ,2.5,3,3.5,4,4.5,5);
$Y=array(0,80 ,160 ,240  ,320 ,360 ,400 ,440,480);
$points = new Matrix(array($X, $Y));
$points->toHTML();
printpoly(polyfit($X, $Y, 6));

echo '<hr />';
?> 